HETEROCLINIC ORBITS, MOBILITY PARAMETERS AND STABILITY 

FOR THIN FILM TYPE EQUATIONS 

R. S. LAUGESEN AND M. C. PUGH 

Abstract. We study numerically the phase space of the evolution equation 

hf — ~{h"'hxxx)x — !3{h™hx)x- 

Here h{x, t) > 0, n > and m e K, and the Bond number B is positive. 

We pursue three goals: to investigate the nonlinear stability of the positive periodic and 
constant steady states; to locate heteroclinic connecting orbits between these steady states 
and the compactly supported 'droplet' steady states; and to determine how these orbits 
change when the 'mobility' exponents n and m are changed. 

For example, when n+l<m<n + 2 we know from the companion article that there 
can be three fundamentally different steady states with the same period and volume. The 
first is a constant steady state that is a local minimum of the energy. The second is a 
positive periodic steady state that is a saddle for the energy and has higher energy than 
the constant steady state. The third is a periodic collection of droplet steady states having 
lower energy than either the positive or constant steady states. Here, we find numerically 
that the constant steady state appears to be asymptotically stable and that perturbing the 
positive periodic steady state in one direction yields a solution that tends to the constant 
steady state while perturbing in the other direction yields a solution that appears to touch 
down in finite time. 

Also, we consider the effect of changing the mobility coefficients, /i" and /i™. We change 
them in such a way that the steady states are unchanged and find evidence that heteroclinic 
orbits between steady states are perturbed but not broken. We also find that when there 
appear to be touch-down singularities, the exponent n affects whether they occur in finite 
or infinite time, ft also can affect whether there is one touch-down or two touch-downs per 
period. 
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1. Introduction 
We study the evolution equation 

(1) ht = -{h^h,.,),-B{h'^h,),, 

where n > and m G M, and the Bond number B > 0. This is a special case in one space 
dimension of the equation ht = —V ■ {f{h)'VAh) — V ■ {g{h)'Vh), which has been used to 
model the dynamics of a thin film of viscous liquid. The air/liquid interface is at height 
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z = h{x,y,t) > and the liquid/solid interface is at 2; = 0. The one dimensional equation 
apphes if the hquid film is uniform in the y direction. 

The fourth order term in equation (P reflects surface-tension-type effects and the second 
order term can reflect gravity, van der Waals interactions, thermocapillary effects or the 
geometry of the solid substrate, for example. Typically 

f{h) ~ /i" and g{h) ~ ±i3/i™ 

as h —>■ 0, where 1 < n < 3 and m G M, ;B > 0, hence our choice of / and g in (|I]) as power 
laws. See the introduction of our companion paper |T^ for references to the modeling and 
mathematical literature. 

In [|r^] we proved linear stability and instability results for the positive periodic steady 
states of (|l]). These results were for zero-mean perturbations that have the same period as 
the steady state, or shorter period — longer perturbations always lead to linear instability 
T8| . We summarize the linear stability and instability results in the bifurcation diagrams 



in the main points are roughly that positive periodic steady states are linearly unstable 
ifm — n<OoTm — n> .80, are linearly stable if < m — n < .75 (albeit with a zero 
eigenvalue arising from the translation invariance of (|l]) ), and could be either stable or 
unstable if .75 < m — n < .8. The case m — n = is degenerate. 

For the linearly stable steady states, it is natural to ask: do numerical simulations show 
asymptotic stability of the steady state (modulo translation in space), meaning that a solu- 
tion that starts from a perturbed steady state will relax back either to the steady state or 
to one of its translates? In § ^.4| we show the answer is "Yes" , for m and n with m — n = .5 
and other values. 

For the steady states hgs that are linearly unstable, two natural questions occur. First, 
is the steady state nonlinearly unstable? We find this seems always to be the case, as 
evidenced by the many numerical simulations in Second, do perturbations of the steady 
state subsequently converge to some other steady state, and can we predict what this long- 
time limiting state will be? This is a difficult question. The three most likely candidates for 
a long-time limit are: a positive periodic steady state, the constant steady state /igs, 01 one 
or more compactly supported 'droplet' steady states. (Droplet steady states might have zero 
or nonzero contact angles, but in ||18|, |l^ we have considered only the zero contact angle 
case.) How can one predict which, if any, of these candidates will be the long-time limit? 

One way to make predictions is by comparing the energy levels of these different kinds of 
steady state, using the energy 
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This energy is known to be dissipated by the evolution, ^£{h(-,t)) < 0, and in our companion 
article [|l9l we proved a number of results establishing relative energy levels of positive 
periodic, constant and droplet steady states. For example, we proved in Theorem 6] 
that ifm — n<Oorm — n>l then the positive periodic steady state/igs has higher energy 
than the constant steady state hss, suggesting that perturbations of hss might converge to hss- 
Indeed, we find numerically in §H for a variety of m and n values that many perturbations of 
hss do subsequently evolve towards the mean. In another unstable case, if —2 < m — n < 
orl<m — ra<2 then [|19|, Theorem 7] shows there is a zero angle droplet steady state 
with length less than the period of hss, with the same area as h^s, and with less energy than 
hss- We find numerically in ^ that many perturbations of hss evolve towards touch-down, 
where the solution goes to zero at a point. Our code stops at that time, so we cannot show 
relaxation to the droplet steady state, but this certainly seems a likely eventual outcome. 

Our theorems and simulations thus allow us to predict at least a couple of likely long-time 
limits. But when there is more than one steady state with lower energy than hss, we have not 
been able to predict to where a given perturbation of hss will relax. And we are so far unable 
to predict the amount of translation that might occur in the long-time limit (for example, 
where will the maximum of the limiting steady state occur?). Impressive results on such 
a translation problem have recently been obtained by [jlll for the Cahn-Hilliard equation 
ht = —hxxxx — ((1 ~ 'ih'^)hx)x on the whole real line. 

In this paper we also ask how the evolution is affected by changes in the coefficient functions 
/i" and h^. In particular, we consider the evolution of a fixed initial condition under the 
equation (P for two different pairs {ni,mi) and (^2,^12) of exponents. One can show that 
the steady states of the two evolution equations are the same if mi — ni = 1112 — ^2. In that 
case we ask: given the same initial data, do the two solutions have the same long-time limit? 
Or, can changing the mobility exponents n and m can change the long-time limit? We also 
ask whether changing the mobility parameters can affect the number and type of finite-time 
singularities. 

Before outlining the paper, we note that throughout the paper we consider only zero-mean 
perturbations. This seems reasonable from a physical standpoint, because such perturbations 
correspond to a disturbance of the fiuid that alters the profile without adding additional 
fiuid. Mathematically it is reasonable because the evolution equation (|l|) preserves volume 
for spatially periodic solutions: j h{x,t)dx = J h{x,0)dx for all time t. Thus zero-mean 
perturbations allow the possibility of relaxation back to the original steady state, while 
nonzero-mean perturbations do not. 



Outline of the paper. This paper reports on numerical simulations that are inspired by, 
and extend, the theoretical results in our linear stability paper fl^ and our companion 



4 



paper ||T9l on energy stability and relative energy levels of steady states. The stability 
results of [|18|, |19| are summarized by the bifurcation diagrams in §^ below, and we will 
remind the reader of the relevant results and their implications when describing our numerical 



simulations. So the reader need not digest the earlier papers |T^, |T^] before reading this one, 
although it would help to have those papers at hand. 

The bifurcation diagrams in Section ^ summarize the known results from |T8| , |T9[ on 
stability of steady states. The section also contains a weakly nonlinear stability analysis. 

Section || presents a detailed numerical study of the evolution equation (|1|), for a number 
of different exponents n and m. We focus especially on initial data close to a steady state. 
Our stability and energy level results from ]TB| , ^ lead to many predictions for the behavior 
of the solution, both short and long time, and these predictions are generally borne out by 
our simulations. Strikingly, the period and area of the initial data, the Bond number B and 
the value m — n seem to be sufficient information to reduce the possible long-time behavior 
to just a couple of options. 

We find evidence for heteroclinic connections between different types of steady state (such 
as: periodic to constant, periodic to droplet, and constant to droplet). One would like to 
know how robust the behaviors observed in our numerical simulations are under changes 
in the coefficients and Bh"^. In particular, what happens when the mobility exponents 
n and m are changed in a way that n — m is unchanged, leaving the steady states of the 
evolution unchanged? In Section we give numerical evidence that such changes perturb, 
but do not break, heteroclinic orbits. 

Another question investigated in Section ^ is how changing m and n affects singularity 
formation. For the equation ht = —{h^hxxx)xi there has been extensive computational work 
studying how the choice of n affects the spatial structure of singularities and whether they 
occur in finite or in infinite time Specifically, simulations suggest there is a critical 
exponent 1 < ric < 2 such that ii n > ric then solutions are positive for all time while 
solutions can touch down in finite time if n < ric- For our equation (|l|), we numerically 



estimate Uc in %.2[ We find that its value depends on the difference m — n. Then in § |5.3| , 
we demonstrate that the mobility can also affect the number of apparent singularities per 
period. 

Section ^ discusses our numerical methods. 

Section |^ sums up the paper, and recalls the 'gradient flow for the energy' interpretation 
of the evolution equation (0). 



Figure 1: Four types of steady state. 



2. Terminology, definition of the energy, and review of steady states 

Terminology. We write Tx for a circle of circumference X > 0. As usual, one identifies 
functions on with functions on M that are X-periodic and calls them even or odd according 
to whether they are even or odd on R. 

Positive periodic steady states are assumed to satisfy the steady state equation classically. 
A droplet steady state hss{x) (see Figure ||) is by definition positive on some interval {a,b) 
and zero elsewhere, with hgs G C^[a,b]] we require /igs to satisfy the steady state equation 
classically on the open interval (a, 6), and to have equal acute contact angles: < h'^s{a) = 
— < oo. (Throughout the paper, if a function has only one independent variable then 
we use ' to denote differentiation with respect to that variable: h'^^ = {hss)x-) 

We say a droplet steady state hss has 'zero contact angle' if = h'^si^^) = — /iss(&), and 
'nonzero contact angle' otherwise. A 'configuration' of droplet steady states is a collection 
of steady droplets whose supports are disjoint. 

Definition of the energy. There is a well-known dissipated energy (or Liapunov function) 
for the evolution equation (^. It is defined for i G H^(Tx) to be 



where if is a function with H"{y) = By^ This energy is strictly dissipated: if h{x,t) is 



hence it cannot distinguish between a steady state and its translates. 

A precise definition of linear stability can be found in Appendix A], and of energy 
stability in [0, §2], but those definitions are not needed to understand this paper. 

Brief review of steady states. Here we quickly review the basic facts about steady states 
needed to appreciate this paper. For more on the steady states and their properties, and for 
justifications of the following remarks, see [0, §2.3] and and the references therein. 

We start with a non-constant positive periodic steady state hss ^ C'^iTx) of the evolution 
(HI). We translate h^g so that its minimum occurs at x = 0. The steady state equation for 
(0) integrates to give hssh'ss + Bk^h'ss = C for some constant C. 



(2) 




a smooth solution of (|l]) then {d/dt)S{h{-,t)) < 0, with equality if and only if /i is a steady 
state (cf. §2.1]). The energy, like the evolution equation, is invariant under translation: 



The constant C (the flux) equals zero, by integrating hll + Bk^ "'h'^^ = Ch^^ over a period. 
Hence the steady state satisfies 



(3) 

here r{y) := By'^~^ and 



q := m — n -\- 1. 



This exponent q determines many properties of the steady state, including (usually) its linear 
stability. 

Integrating, we find the steady states have a nonlinear oscillator formulation: 



(4) 



, Bhl - D 

/i„„ H = 

q 



for some constant D, when q ^ 0. For g = the analogous equation is h"g + B\oghss — D = 0. 
This oscillator equation contains three constants: q, B, and D, but B and D can be removed 



by rescaling h^s and x (see |T9|, eq. (7)]), leading to the 'canonical' steady state equation 

ki -1 



(5) 
(6) 



k" + 



0, g^O, 



k" + \ogk = 0, q = 0. 



Every positive periodic steady state hss can be rescaled to such a function k = ka with 
k'^{0) = 0, where we write a = k{0) G (0, 1) for the initial value. Conversely, for each g G M 
and a G (0, 1) there exists a unique smooth positive periodic fc„ satisfying equations (||-§|) 
and with ka{0) = a, k'^{0) = 0. The same holds for a = when q > —1, although ko may 
be only C^-smooth at x = 0, where /co(0) = (see Theorem 3.2]). 

We write 



P = P. 



P{a) 



and 



A = Ar 



A{a) 



respectively for the least period of ka and for the area under its graph, A = ka{x)dx. 



Note that Pa and Aa approach 27r as a — > 1. As seen in p!8| , |T9|, the function 

E{a) := P{af-''A{ay^^ = P{a)^[A{a)/P{a)]''-^ 
plays a large role in the stability theory of the steady states, in part due to its rescaling 



mvariance: 



(7) BPt'At' = P{af-'A{ar\= E{a)), 

where Pss and Ass are the period and area of the steady state hss- 
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3. Bifurcation diagrams and weakly nonlinear analysis 



m 



In this section we present bifurcation diagrams that encode the hnear stabihty information 
ll, . Then we relate our bifurcation diagrams to a weakly nonlinear analysis near the 



constant steady state. 

3.1. Bifurcation diagrams. Figure ^ gives bifurcation diagrams for various values of q. 
These g-values are representative members of the intervals (— oo, 1), (1, 1.75), (1.75, 1.794), 
and (1.795, oo). (The value 1.794 approximates a critical exponent; see §5.1] for further 
details.) 

In all these diagrams we are considering the 27r-periodic problem with B = 1. We construct 
the diagrams as follows. Given q, we first compute the rescaled steady states ka for a range 
of a G (0, 1). For each ka, we use its period Pa to determine a constant D in the rescaling 
il9| , eq. (7)] that yields a positive steady state hss of period Pss = 27i. We then plot the 
amplitude of hss versus the scale invariant quantity E = BPj^^'^A'^^^, determining the linear 
stability of the steady state by the results in [T^, §3.2], particularly [|18], Theorem 9]. The 
conclusions on linear stability summarized in the figure are rigorously proved in except 
for 1 < g < 2, in which case we are relying on a combination of analytical and numerical 
results. 



£/ = -3 




Cf = .5 



£J= 1.75 



q= 1.775 



q=1 



q= 1.76 



q = 2.5 



Fig. 2: The x-axis is BP^^ '^A^^ ^ and y-axis is [hmax — 
unstable; dotted: linearly neutrally stable; solid: linearly stable. 



i)/2. Dashed: linearly 



The horizontal axes of these diagrams also show the linear stability of the constant steady 
state, with respect to zero-mean perturbations, where the constant steady state is consid- 



ered as having period i^s and area Ass- This linear stability information is taken from ||19 
Theorem 10]. 
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Qualitatively, the diagrams change continuously as q increases. Choosing a period other 
than 27r or a Bond number other than B = 1 simply dilates the y-axis of the diagrams. We 
discuss the diagrams further when we present numerical simulations in §|. 

3.2. Weakly nonlinear analysis. In the following, we sketch the weakly nonlinear analysis 
for the evolution equation (|l|). This was done for q = —3 in |26, §2.3]. 



In short, we consider physical parameters such that the constant steady state has one 
mode which is barely linearly stable or is barely linearly unstable, while all other modes are 
strongly damped. This yields a separation of timescales which allows one to find a reduced 
representation of the PDE in terms of an ODE governing the amplitude of the unstable 



mode. We refer the reader to §5.1 of |20 



Let h be an X-periodic solution of ht = —{h"'hxxx)x — B{h"^hx)x with mean value h. We 
rescale the solution to have period 2tt and mean value 1, and we also rescale time: 



27CX — , , , /27r^^ 



C = hr] = h, and t' = j t. 

The rescaled evolution equation is rjf = — o"('7"^7ccc)c ~ /^(^™^c)c where 

ZTT 



Linearizing about h = 1, we find that linearly unstable modes exist if < k < k^. Proceeding 
in the usual manner, we introduce a small parameter 6 that corresponds to moving through 
the critical wave number kc = 1: 

^ = 1 + Q5' 
a 

where Q = ±1. We then introduce a slow time-scale r = and expand the solution in 
orders of 6: r) = 1 + 5?7i(C, t) + 6'^ri2{C, r) + 6^ri3{(, r) + 0(5^). For simplicity, we assume 
the solution is even. By the usual arguments, ?7i(C7 t) = A{t) cos(C), i]2i(, t) = B{t) cos(2^), 
and rj^iCr) = C(r) cos(3C). Putting this ansatz into the evolution equation and expanding 
in orders of 6, we find there are no 0(1) or 0{6) terms. At 0(5^) and 0{6^), one determines 
the amplitudes B{t) and C{t) in terms of A{t) which, in turn, satisfies 

^ = QaA{T) - nAirf where k = -(g - 1)(1.75 - g). 
dr 6 

The dynamics of the amplitude A{t) depend on the signs of the Landau constant k, and 
of the linear term. If k > then for Q = 1 the constant A{t) = is linearly unstable 
and A{t) saturates to the linearly stable amplitude Ac{t) = o j k. This corresponds to 
a supercritical bifurcation. If k < then for Q = —\ the constant Air) = is linearly 
stable and the steady amplitude Acij) = —o is linearly unstable. This corresponds to 
a subcritical bifurcation. Since cr > 0, 

1 < g < 1.75 =^ supercritical bifurcation and g < 1, 1.75 < g =^ subcritical bifurcation. 
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The Landau constant k is determined by the mean value h, the period X, and the ex- 
ponents m and n, but not the Bond number B. Subcritical bifurcations are often seen in 
systems that can have finite-time pinching (rupture) singularities e.g. §3.2], |^ §IV]. 

The above weakly nonlinear analysis is consistent with our analytical and computational 
results. Specifically, recalling the bifurcation diagrams in Figure ^ we note that positive 
periodic steady states are linearly unstable for g < 1 and q > 2. For 1.75 < q < 2, when 
there is only one positive periodic steady state it is linearly unstable. When there are two, 
the one of smaller amplitude {hmax — hmin)/2 is linearly unstable. For 1 < g < 1.75 the 
positive periodic steady state is linearly stable. 



4. Simulations, and heteroclinic orbits connecting steady states 

We now numerically simulate solutions of the evolution equation ([^), for a wide range 
of initial data near steady states. This has not been done before. The solutions obtained 



display a great variety of stability and long-time behaviors. Our stability theorems ||18| , [19 
often allow us to predict the numerically observed short-time behavior, and our theorems 
on the energy levels of steady states |19| often allow us to guess the long-time limit of the 



evolution. As part of this we predict (and find strong evidence for) heteroclinic connections 
between certain steady states. 

We expect our numerical investigations of the power law evolution (|1]) will provide re- 
sources, ideas and motivation for researchers studying ht = —{f{h)hxxx)x — {g{h)hx)x with 
non-power law coefficient functions / and g. There are some such numerical studies already. 



For example, the papers [|T^, ^ consider an / that is degenerate (/(O) = 0) and g^s that 
are not power laws, and there is a large literature on the Cahn-Hilliard equation (for which 
/ = 1 and (? is a quadratic). 

Recall the steady states depend on the parameter 

q = m — n + 1. 

We take seven values of q: 

q = -3, 0.5, 1, 1.5, 1.768, 2.5, 4, 

representatives of the intervals {(— cxo, — 1], (— 1, 1), (1, 1.75], (1.75, 1.79), (1.8, 3), [3, oo)} in 
which our theorems suggest the solutions will display distinct behaviors. For the q = —3 
case, we take n = 3 and m = — 1, making (|l]) a 'van der Waals' equation previously studied 
by other authors. Otherwise, we take n = 1 and m = q. In any event, we find in §^ that our 
numerical simulations are not greatly affected if we change n > and m in a manner that 
keeps q fixed {i.e. that keeps m — n fixed). 
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4.1. q = —3: the van der Waals case. 

Characteristic features for q G {—oo, —1]: positive periodic steady states are linearly un- 
stable, and there are no droplet steady states with acute contact angles. (See bifurcation 
diagram l^a and §2.2].) 



4.1.1. q = —3. Perturbing the positive periodic steady state. First we explain how to find a 
steady state having period P^s = 27i and having some specified value for the area (or volume) 



Af.f. — /q hss{x) dx. 

Given the Bond number B (a physical parameter), a positive steady state satisfies (^: 

Ks + — = 

q 

for some constant D. The problem is to find the value of D for which there is a steady state 
of period 2-7? and area A^s- For q < —1, if A^s > 2tc B^^^^~''^ then there exists a non-constant 



positive periodic steady state h^s with period 27r and area Ass (see [0, §5.1]). The rescaling 
|19| , eq. (7)] then implies there is an admissible value of D. 

For simplicity, instead of starting with the Bond number B and finding this admissible 
values of -D, we instead fix D = 1 and determine the interval of admissible B values. That is, 
we choose a G (0, 1) and consider the steady state ka of the rescaled equation, as in §||. Its 
period P{<y) together with D = 1 and Pss = 27r then uniquely determine B by |l^, eq. (24)], 
and hence the steady state hss by |jl9|, eq. (7)]. (Choosing a different a would yield a different 
B and hss-) Finally, we locate a nearby 'finite-difference steady state' on N meshpoints (see 
§ |67^ ) and study numerically its stability under perturbation. 

For q = —3 we carry out this construction with a = 0.2145 and period Pa = 16.32, 
resulting in Bond number B = 0.003259, and a non-constant positive periodic finite-difference 
steady state hss of least period 2tt and with area A = 6.884. Note the positive periodic steady 
state is linearly unstable, by bifurcation diagram ^ with BPs!^~''Aj~^ = .08930. 

Remark: For the rest of whenever we refer to a 'positive periodic steady state', we 
implicitly mean a 27r-periodic finite-difference steady state that has its minimum at a; = 
(see § |6.3| )- Also, we always find our periodic steady states as above, by fixing D = 1 then 
choosing a G (0, 1) and then determining B and Ass- The only exception is for the q = 1.768 
simulations discussed in § |45| . 

Now that we have a steady state for q = —3, we use a perturbation of it as initial data in 
the equation 

(8) ht = -{h^h^^^)^ - B{h-^h^)^. 

Here n = 3, m = — 1 and q = m — n + l = —3. We study this equation for the rest of § [4.1| . 
It was proposed by Williams and Davis [0] to model a thin liquid film with net repulsive 
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van der Waals interactions, and more recently it has been studied by Zhang and Lister ^7 
and by Witelski and Bernoff |2^, ^ . 



• Even perturbations. 

First we perturb hss with the even zero-mean perturbation ^eh'^^ and numerically study 
the resulting solution. 

The steady state /igs is linearly unstable and since the perturbation +£h"s lowers the 
maximum of hss and raises the minimum, one might hope the resulting solution would 
converge to the constant steady state h = hss- If this happens for all small e, then this 
would be strong evidence for existence of a heteroclinic orbit connecting hss to the constant 
steady state. There are a number of theoretical reasons to suspect such heteroclinic orbits 
exist: (i) hss is energy unstable in the directions ±/iss by [0, Theorem 2], (ii) the energy of 
hss is higher than that of the constant steady state hss by Theorem 6] (also observed 
numerically by Witelski and Bernoff |^6|, §3]), and (iii) the constant steady state is a local 
minimum of the energy S by [|19|, Theorem 10]. [To deduce (iii) from |[19| , Theorem 10] 
requires Bhss < 47r2, which we now establish: Bhss = BX^-^A^-^ = E{a) for 

some a e (0, 1), by rescaling as in (|^), and E{1) = Att^. So we want to show E{a) < E{1)] 
this holds because E' > by |]18|, Theorem 11], since g < 1.] 

To seek evidence for a heteroclinic connection, we start with initial data hss + lO^^/igg. 
Here, we've normalized h'ss to have L°° norm 1; all perturbations are similarly normalized. 
We have also considered e smaller than 10~^, in most of the simulations below; we found 
that these smaller perturbations resulted in qualitatively the same behavior, until one hits 
the level of roundoff error. 




Fig. 3: q = — 3, n = 3. Dashed: initial data hs, 
and solution relaxes to the mean. 



10 /igg. Local extrema stay fixed in space 



The steady state hss has very large curvature at its local minima, and so we need a large 
number of meshpoints to resolve the initial data hss + 10~^/iss with spectral accuracy. We 
find that for a solution on [0, 27r) we need 2048 meshpoints: hss{x) = X]fc=-io23 '^fc ^^P(^^^) 
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has amplitudes that decay to the level of round-off error (a^ ~ 10 for k ~ 1024). Figure ^ 
shows the evolution of the solution; the solution relaxes to the constant steady state. (This 



was shown previously in Figure 4b].) We use the adaptive timestepping described in § |6.2| . 
At first sight, this would seem unnecessary since the solution is becoming more regular as it 
evolves. However, there is a short transient during which hmm{t) decreases. The timestep 
initially decreases to accurately track this fast behavior and then increases as the solution 
relaxes to the mean. We use adaptive time-stepping throughout our work since we almost 
always observe such a short transient. Also, in a number of cases, the solutions become less 
smooth (curvatures increase) as time passes, requiring refinement later in time. 

In the opposite direction, the perturbation —eh'^^ raises the maximum and lowers the 
minimum of the steady state h^s. Since h^g is energy unstable, we might expect the solution 
to subsequently converge to a droplet steady state or to a configuration of droplet steady 
states. From fl^, §2.2], if such a droplet exists it must have 90° contact angles, though we 



have not discussed such steady states here in this paper or in ||17| , [18 



Our numerical simulation of the solution with initial data hss — 10~^/igg shows that, after a 
short transient, the minimum height of the solution decreases in time, appearing to decrease 
to zero in finite time. As the minimum height decreases, the curvature increases, requiring 
that after some time the number of meshpoints be increased to keep the solution spectrally 
resolved. We do this as follows. We compute the solution with 2048 meshpoints until the 
computation stops {hmmif) = 0, see §|6l4|.) We look at the power spectrum of the solution 
and choose a time right before the active part of the power-spectrum is reaching the Nyquist 
frequency (see right plot of Figure |^). That is, we find the last time at which the 1023rd 
Fourier amplitude of the solution is at the level of round-off. We take the solution at this 
time and compute its Fourier coefficients, defined for wave numbers — n/2 + 1<A;<n/2 — 1 
where N = 2048. We pad by zeros, extending the Fourier coefficients to be defined for 
wave numbers — N + l</c<N — 1, and then compute the inverse Fourier transform. This 
yields a function on 2n meshpoints that is indistinguishable from the solution at that time, 
to the level of round-off. Using this function as initial data, we continue the computation 
on 2n meshpoints, repeating this point-doubling process whenever the solution becomes 
unresolved. 

In this way, we computed a resolved solution to time t = .0472496406249, the time when 
the 32,768 meshpoint solution became unresolved. The top left plot in Figure]^ presents the 
evolution of the solution near x = 0. As before, the local extrema are fixed in space, with the 



solution appearing to touch down at one point per period. (This was shown previously in p6 



Figure 4c].) We did not design the code to study the formation of finite-time singularities; 
the solution has decreased by only a factor of 6.04 due to limitations of our uniform-mesh 



code (see § |6.4|) . Computing the derivative of the solution, we find that its maximum 
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and minimum values grow as time passes, as in the bottom left plot of Figure ^. These 
extremum points of move in time, appearing to converge to s = as the singular time 
approaches. This is consistent with a solution that touches down with 90° contact angles in 
finite time. The right plots in Figure ^ show the final resolved solution, which suggests 90° 
contact angles are developing. 
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Fig. 4: Left top: dashed line, initial data /igs — ^^~^h'^s'j solid lines, solution h at later times. 
Local minimum is fixed in space and decreases monotonically after short transient. Left bottom: 
dashed line, initial slope; solid lines, hx at same times as above (||/ia;|loo increases in time). Right 
top: solution at final resolved time, t = .0472496406249. Right bottom: close-up near x — Q. 



The work of Zhang and Lister §5] on similarity solutions suggests that t) ~ 
B^^^{tc — ty^^H(x/{tc — tY^^) as touchdown approaches; here tc is the time of touchdown and 
if is a particular positive function with H(r]) ~ (0.807) for large r]. Our computations 
are consistent with the above ansatz. Further, if we make the ansatz then we can estimate 
tc, since taking the ratio of the computed values of h{0,t) at two late times ti and ^2 gives 
the value of (tc — tiY^^ /{tc — ^2)^^^, from which tc can be determined. We find tc is slightly 
larger than the final resolved time. See p5|, ^ for more on the similarity solutions to (H). 



The most plausible future behavior for the solution in Figure ^ is that the solution might 
touch down in finite time and become a nonnegative weak solution. Then it might relax, 
as a weak solution, to a droplet steady state. Alternatively, since the energy of the original 
periodic steady state is higher than the energy of the constant steady state, the solution might 
perhaps touch down in finite time and then at some later time become positive and smooth 
again, ultimately relaxing to the constant steady state. This is certainly possible since the 
solution shown in the right of Figure ^ has higher energy than the constant steady state. 
However, we did not write our code to study finite-time singularities or weak solutions, and 
consequently we cannot distinguish which (if any) of the above options might be happening. 

We close with a graphic demonstration of the kind of spurious effects that can occur if the 
solution is not spectrally resolved. In the left plot of Figure H we present the 8192 meshpoint 
solution that starts from the same initial data as in Figure The right plot of Figure |^ shows 
the corresponding power spectra: as the solution evolves, higher and higher frequencies are 
needed to resolve the ever-sharpening local minimum, until at the final time resolution has 
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been lost. At this time the solution has multiple oscillations. We found those oscillations 
then grew and the solution seemed to touch down in finite time with two droplets per period, 
one large and one small. But the small droplet is a numerical artifact, in view of the absence 
of small droplets in the resolved solution shown in Figure H. 
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Fig. 5: Left: evolution near x — oi the solution with 8192 meshpoints. Right: the cor- 
responding power spectra. The power spectrum reaches the Nyquist frequency soon after 
t = .04716, and the solution loses resolution. The profile on the left with multiple local 
minima is after t = .04716. The profiles with a single local minimum are spectrally resolved. 

• Odd perturbations. 

We now turn to odd perturbations. An obvious choice is the perturbation eh'^g, it arises 
naturally from a leftwards translation of hss, since hss{x + e) ~ hss{x) + eh'g^{x). This 
perturbation shifts the local minimum of the initial data from a; = to a point slightly to 
the left of X = 0, and it also lowers the minimum value since for small x < 0, 

/iss(x) + eh^x) = /iss(O) + eKM^ + 0{x^) < h,{0). 

Given initial data hgs + lO"^/^^^, we expect the solution either to approach zero at some 
point or to relax to the constant steady state, since hgs is energy unstable in the direction 
±/igg by Theorem 2]. The simulation confirms this expectation, yielding a solution that 
touches down in a fashion very similar to the solution shown in Figure ^. One difference is 
that the location of the local minimum is no longer fixed in space; it moves as hmin{t) —>■ 0. 
The similarities are that the minimum value of the initial data is very close to zero, relative 
to the bulk of the initial data, and that the local minimum of the solution appears to decrease 
to zero in finite time while the bulk of the solution is essentially unchanged from the initial 
data. 

• Random perturbations. 

The equation (|^) is translation invariant so there is nothing to distinguish any particular 
point in space. Up until now, we used either even or odd perturbations, making x = appear 
somewhat distinguished. To check the degree to which the behaviors described above depend 
on the choice of perturbation, we performed a number of runs with random perturbations. 

The random perturbations are constructed as follows. The 2048 meshpoints can resolve 
1023 frequencies, so we chose {ak} and {(pk} to be two sets of 1023 uniformly distributed 
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random numbers in [0, 1] and defined the perturbation 



1023 

(f){x) = Ofc exp(— .036A;) cos(/cx + 2710^ j 

k=l 

This perturbation has zero mean, and has random amphtudes and phases at each wave- 
number. The decay rate 0.036 is chosen so that the amphtudes at wave numbers k = 900 
and higher are at the level of round-off error, so that the initial data is spectrally resolved. 
We then divided by its norm. 

We found that all the solutions resulting from applying such a random perturbation to hss 
either relaxed to the constant steady state or else appeared to touch down in finite time. An 
intuitive rule might be that if the minimum value of the initial data is greater than that of 
then the solution should relax to the constant steady state and if the minimum value is 
less than that of hss then the solution should touch down in finite time. In practice, we found 
that most solutions respected this intuition, and the solutions that appeared to touch down 
in finite time had their gross dynamics as described earlier (their finer dynamics concern the 
position of the local minimum as a function of time). However, there were exceptions — 
hardly surprising since the evolution equation is a fourth order PDE, not second order, and 
the intuitive rule has the flavor of a comparison principle. 

4.1.2. q = —3. Perturbing the constant steady state. Consider the constant steady state 
h = h, for some Bond number B. If Bh'' ^ < 1 then /i is a strict local minimum of the 
energy by Theorem 10], and is nonlinearly stable. In this case a non-constant positive 
periodic steady state hss exists with period 27r and mean value h by |jl8|. Theorem 12], since 



< Bh'^ ^(27r)^ < (27r)^; it is linearly unstable (see bifurcation diagram ^ja). If Bh'^ ^ > 1 
then the constant steady state h is a. saddle point for the energy and we expect it to be 
unstable. In this case there is no stable positive periodic steady state to which a perturbation 
could converge (see bifurcation diagram ^), suggesting that the solution will touch down in 
finite or infinite time. 
• Cosine perturbations. 

We take B = 1. Consider the constant steady state h = 2, a. local minimum of the energy 
since Bh"^ ^ < 1. Here the constant steady state is linearly stable and there exists a linearly 
unstable positive periodic steady state, by above. Given initial data 2 — 10"'^ cos x, the 
solution appears numerically to relax to the constant steady state. 

On the other hand, the constant steady state /i = 1/2 is linearly unstable and is a saddle 
point of the energy. Indeed, for initial data 1/2 — 10^'^ cos x we find a solution that appears 
to touch down in finite time. Figure ^ shows this evolution over two periods. The top plot 
shows the short-time dynamics: the local minimum decreases, while the local maximum 
increases for a while. The top then "flattens" and two local maxima form one to each side 
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Fig. 6: q — — 3,n = 3. Dashed: initial data 1/2 — 10 '*cos(x). The local minimum is at x = 
at all times and decreases monotonically to zero. 

of the flat region. Tlie bottom plot shows the later-time dynamics: the solution appears 
to touch down at one point per period and continues to have two local maxima per period. 
The final profile presented does not look like any known steady state, and so we expect the 
solution to continue evolving after touching down. 
• Random perturbations. 

Using random perturbations, we verified that the above results are robust: the steady 
state h = 2 is asymptotically stable and = 1/2 is unstable. 

Incidentally, we verified that the evolution is exponential in time near the periodic and 
constant steady states. This is consistent with the nonlinear behavior being dominated by 
the linear theory when the solution is sufficiently near a steady state. We found there was 
a short transient before the exponential behavior began, suggesting that the direction h'^^ 
is near but not equal to the first eigendirection. During the transient time, the solution is 
locating this eigendirection. 

4.2. q = 0.5. 

Characteristic features of q ^ (—1, 1).- positive periodic steady states are linearly unstable. 
A 'Mountain pass' scenario can occur — the energy of the non-constant positive periodic 
steady state is higher than the energies of the constant steady state and a zero- contact angle 
droplet steady state. (See bifurcation diagram and remarks after Theorem 11].) 

We take n = 1, m = .5, g = 0.5 and compute solutions of 



4.2.1. q = 0.5, perturbing the positive periodic steady state. We rescale a steady state 
with minimum height a = 0.03000 and period P{a) = 6.049. This yields a Bond number 
B = 0.9628 and positive periodic steady state hss with period 27r and area Agg = 7.980. Note 
hss is linearly unstable, by bifurcation diagram ||b with BPgl~'^A1~^ = 35.59. 
• Even perturbations. 



■t — ~ 
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As in § |4.1.1| , we perturb hss with ±10^^h'^^. 

For the initial data hss + 10~^/;.ss, we expect to see evidence of a heterochnic connection 
to the constant steady state hss, for the same reasons given for the q = —3 case, and indeed 
our simulation turned out to be very similar to that shown in Figure |^. 

Next, for the initial data hss ~ 10~^/igg we find the solution appears to touch down in 
finite time. Like the q = —3 simulation, hmin{t) is located at x = and, except for a short 
transient, decreases monotonically in time. Figure |^a presents h at the final resolved time. 





Fig. 7: q — .5, n — 1. Left: h at final resolved time, t — 7.762540625, shown over two 
periods with 8192 meshpoints per period. Right top: close-up of h near pinch-ofT; minimum 
decreases monotonically after a short transient; dashed line is initial data. Right bottom: 
close-up of hx at same times. 



Since — 1 < g < 1, [0, Theorem 7] tells us there exists a zero-angle droplet steady state hss 
that has the same area as hss, has length less than 2it, and has lower energy than hg^- And 
indeed Figure |^a shows a profile that appears to be close to a zero-contact-angle solution. As 
further evidence of this, in the top plot of Figure 0b we present a close-up of the evolution 
near the touch-down point. In the bottom plot we present h^ at those times: the slope does 
not appear to be forming a jump discontinuity, which it would have to be doing were the 
solution converging to a solution with nonzero contact angle. 

• Odd perturbations. 

The odd perturbation 10~^/igg yields an evolution qualitatively like that of g = —3 above. 
The solution appears to touch down in finite time, with the late-time behavior much like 
that shown in the right top plot of Figure |^. 

• Random perturbations. 

For random perturbations, we observed the same type of dynamics as seen in the q = —3 
case. Some perturbations led to solutions that relaxed to the constant steady state, and 
others yielded solutions that appeared to touch down in finite time with a late-time evolution 
like in the right top plot of Figure ^ 

Mountain pass: Our simulations above numerically confirm for q = .5 the following 'mountain 
pass' scenario, which is possible whenever — 1 < g < 1. 
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Choose a positive periodic steady state hss such that B{27t)^^'^ A^^^ lies between EQ{q) and 
477^ (see Figure 6] and refer to the existence resuh |jl8|, Theorem 12]). This positive 
periodic steady state is hnearly unstable and has higher energy than the constant steady 
state h = hss, by 0, Theorem 7] and [|19], Theorem 6]. It also has higher energy than the 
zero-angle droplet steady state hss by ||19|, Theorem 7]. Further, the constant steady state 
hss is linearly stable and is a local minimum of the energy, by |T9| , Theorem 10]. 

Thus the positive periodic steady state hss appears to sit at a 'mountain pass' between 
the constant steady state hss (which lies at the bottom of an energy well) and the droplet 
state hss- If in addition B{2Tr)^~'^A1~^ > L{q) (see [|19|, Figure 6] for L{q) ) then the constant 
steady state has higher energy than the droplet, by []T^ Theorem 11]. But regardless of 
that, one would expect that perturbing hss in one direction would lead to a solution that 
converges to the constant steady state while perturbing in the other direction would lead to 
a single droplet. Our numerics for g = .5 are all consistent with this expectation. 




Fig. 8: q ^ .5,n ~ 1. Dashed: initial data 1/2 — .0001 cosx. The local minimum and maximum 
are fixed in space and, after a short transient, decrease (increase) monotonically. 

4.2.2. q = 0.5, perturbing the constant steady state. As in the q = —3 case, we take B = 1 and 
consider perturbing the constant steady states h = 2 {Bh'^ ^ < 1) and h = 1/2 {Bh^ ^ > 1). 
The first steady state is linearly stable while the second is linearly unstable. We perturbed 
the two states with a range of perturbations, and found that all perturbations oi h = 2 
relaxed to /i = 2, while all perturbations of /i = 1/2 led to apparent finite-time touch-down. 

In Figure ^ for example, we plot the evolution of the solution with initial data 1/2 — 
.0001 cos(a;). The extrema are located at x = and 27r and, after a short transient, increase 
(decrease) monotonically. The evolution shown here is very standard; we observed this type 
of behavior more often than the type shown in Figure ^ 

4.3. q= 1. 

Characteristic features for q = 1: all positive periodic steady states are linearly neutrally 
stable. (See bifurcation diagram |c and [|l^. Lemma 4].) 
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Here we consider the q = 1 case, for which the non-constant positive periodic steady states 
are neutrally stable. We take n = m = 1 and compute solutions of 



Goldstein et al. [0, Fig. 3a] found that fairly large multi-modal perturbations of such steady 
states relax to steady states. They found the solution may relax to a different steady state 
than the one of which it was initially a perturbation. 

For B = 2(1 — cos(Aa;))/(Aa;)^ = 1 — 0{Ax^), one finds a finite-difference steady state 
by sampling a + 6cos(x) + csin(x) on a Ax-uniform mesh (see §6.3| ). In the left plot of 
Figure we present two simulations confirming that each positive periodic steady state is 
nonlinearly stable and that a small perturbation of a positive periodic steady state con- 
verges to a (potentially different) positive periodic steady state. In the top left plot, we 
present the evolution from initial data 1 — .8cos(a;) + .3v{x) where f is a zero-mean ran- 
dom perturbation. In the bottom left plot, we present the evolution from initial data 
1 - .8cos(x) - .19exp(-100sin2(x/2)) + .19 exp(-100 sin2((x - vr)/2)). In both cases, the 
solution relaxes to a positive periodic steady state with an amplitude close to .8 and a local 
minimum close to x = 0. We find that the smaller the perturbation, the closer the long-time 
limit is to the original steady state, numerically demonstrating nonlinear stability. We have 
no rule for predicting the amplitude of the long-time limit and, unless the perturbation is 
even, we have no way of predicting the position of the local minimum. 

Since these simulations suggest that the non-constant positive periodic steady states are 
nonlinearly stable, one might guess that one cannot find a solution that touches down in 



finite time. This is certainly what we observed for q = 1.5 and q = 1.768 in § [4.4| - [4.5| . And 
as the bottom left plot of Figure |^ suggests, initial data that has a sharp local minimum will 
likely not evolve towards touch-down; the local minimum will retract in time, as expected 
for a solution of a surface-tension driven flow. But initial data that is very fiat near its local 
minimum, such as ho{x) = .7 — .8cos(a;) + .19cos(2a:;), does appear to lead to touch-down 
in finite time, as shown in the top right plot of Figure ^. The bottom right plot shows the 
local evolution near the touch-down point. 

The q = 1 case remains mysterious in many ways, because there are infinitely many 2n- 
periodic steady states all having the same mean value; in the q ^ 1 case there are at most 
two. 



Note: numerical simulations for q = 1 were earlier presented with m = n = 1 m |T^, with 



m = n = 2 in [ll5| , §8], and with m = n = 3 in [g]. The latter two articles do not consider 



Bond numbers for which periodic steady states might be observed. 
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Fig. 9: q — \,n—\\ dashed line: initial data. Top left: hQ{x) — 1 — .8cos(a;) + .3w(a;) where 
I) is a random perturbation. Bottom left: ho{x) = 1 — .8cos(a;) — .19 exp(— 100 sin^(x/2)) + 
.19exp(— 100sin^((a; — 7r)/2)). In both cases, solution relaxes to a steady state. Right: 
ho{x) = .7 — .8 cos(a;) + .19 cos(2a::) with minimum value .09. Extrema are fixed in space and 
change monotonically (after a short transient), with /imm(0 touching down in finite time. 
The final resolved solution is on 32,768 meshpoints with hmin — 1-02 10^^. 

4.4. q = 1.5. 

Characteristic features for q G (1, 1.75].' positive periodic steady states are linearly stable. 
(See bifurcation diagrams ||d-e.) 

We take n = 1, m = 1.5, q = 1.5, and compute solutions of 

hf = —{h^hxxx)x ~ ^{h^'^hx)x- 

4.4.1. q = 1.5. Perturbing the positive periodic steady state. We rescale a steady state 
with minimum height a = 0.2145, period Pa = 6.453. With D = 1 this yields a Bond number 
B = 1.083 and positive periodic steady state hss with period Pss = 2tt and area Ass = 5.516, 
and with minimum at a; = 0. This positive periodic steady state is linearly stable, unlike for 
the other g-values considered so far; see bifurcation diagram ^d with BP,l~''Aj~^ = 40.07. 
• Even perturbations. 

We perturb hss with ±10~'^h"^, and expect to see the solution relax back to hss, since hss 
is linearly stable and since translation in space is ruled out by the evenness of hss and of 
our perturbation. Numerically, we indeed observe relaxation back to hss'- the extrema are 
at a; = and n and, after a short transient, they move monotonically to the maximum and 
minimum values of hss- 
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Fig. 10: q — 1.5, n — 1. Dashed line: initial data /igs + -S/iss- Isolated solid line: h^^. Solution 
relaxes to a translate of hss] extrema move in time. 
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• Odd perturbations. 

We now perturb /igs with .5h'^^, which is a large perturbation relative to those discussed 
above. In view of the linear stability and the fact that every perturbation increases the 
energy, by Theorem 5], we expect to observe relaxation back to a translate of hss- (Note 
the energy and our linear stability results are insensitive to translation.) Convergence back to 
hss itself (with no translation) seems unlikely since the odd perturbation breaks the evenness 
of the initial data. 

Figure |T0| displays the evolution of the solution with initial data hss + -S/iss! ^he solution 
certainly seems to converge to a translate of hss- We verified this type of behavior for a range 
of g G (1,1.75]. 

• Random perturbations. 

Using random perturbations, we verified the robustness of the above results: the periodic 
steady state hss is asymptotically stable, up to translation. 




Fig. 11: q — 1.5, n — 1. Top plot: dashed line is initial data hss — .0001 cos(a;/2); heavy solid 
line is final resolved solution a.t t — 61.7305, with 4096 meshpoints. Bottom plot: close-up of 
final resolved solution. 



4.4.2. q = 1.5. Perturbing the positive periodic steady state with longer perturbations. We 
have demonstrated above that the positive periodic steady state hss is asymptotically stable 
to small perturbations of the same period, 2tt. However, it is linearly unstable to zero-mean 
perturbations of longer period — An, Git and so on — by ||1^, Theorem 1]. 

For initial data hss ~ .0001 cos(a;/2), i.e. a perturbation with period 47r, the top plot of 



Figure |TT] presents the evolution of the solution. The steady state hss has minimum heights 
at X = 0,27r, while the perturbation — .0001 cos(x/2) is even about x = In and decreases 
the initial value at x = and increases the value at x = 27r. The solution appears to touch 
down in finite time, though it does not do so at x = 0; also the solution does not appear to 
be converging to a single droplet. The bottom plot of Figure shows a close-up of the final 
resolved solution. The smaller droplet is not close to a steady droplet, since it contains a 
local minimum within itself — an impossibility for a steady droplet. We expect the solution 
would continue to evolve as a nonnegative weak solution, relaxing either to a single steady 
droplet or to some (unknown) configuration of steady droplets. 
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We also considered a number of random 4'7r-periodic perturbations, and always found 
that the solution appears to touch down in finite time with one large droplet flanked as in 
Figure |ll] by a smaller profile which contains a local minimum within itself. 

4.4.3. q = 1.5. Perturbing the constant steady state. We take i3 = 1. As always, the constant 
steady state hss = h is a strict local minimum of the energy if Bh'^ ^ < 1 (or /i < 1), and is 
a saddle point if Bh'' ^ > 1 (or /i > 1). 
• Cosine perturbations. 

We consider h = .5, 1.05, and 2. We expect h = .5 will be stable to all perturbations, and 
that h = 1.05 and h = 2 will be unstable to some perturbations. 

For h = 1.05, there exists a linearly stable positive 27r-periodic steady state with mean 
value 1.05, by bifurcation diagram |d, using that E := Bh'^ ^(2vr)^ G (An'^ , Eq^q)) when 
B = l,q = 1.5 and h = 1.05, since Ef){q) = 40.67. (See ||17| , §3.1.2] for the formula for 
Eo{q).) A perturbation of h = 1.05 might converge to this positive periodic steady state, 
especially since: hgs is linearly stable; hss should have lower energy than h, as discussed after 
||T9| , Theorem 6] ; and there is no zero-angle droplet steady state with length less than 27r and 
the same area as h, by |19|, Theorem 8]. (Note: we cannot exclude that the solution might 
converge to a droplet steady state with nonzero contact angle.) 

On the other hand, for h = 2 there is no positive steady state with least period 27i and 
mean value 2, since Bh'^ ^(Svr)^ ^ (47r^, £'o(q')). (Or see bifurcation diagram ^.) For this 
reason, we expect that perturbations of h = 2 should lead to solutions that converge to a 
configuration of steady droplets. 




Fig. 12: q = 1.5, n = 1. Dashed line: ho{x) = 2 — .0001 cos(a;). The maximum points are fixed; 
the minimum is initiaUy at a: = but then sphts into two minima which move apart in time as 
they decrease. 

We take initial data h — 10~^cosx. We find that if h = .5 then the solution relaxes to 
the constant steady state and if h = 1.05 then the solution relaxes to the positive periodic 
steady state. In both cases, the dynamics are very natural, somewhat like Figure ^ (but 
run in reverse when h = 1.05), with the local extremum points fixed in space and the 
corresponding extremal values evolving monotonically in time, after a short transient. For 
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h = 2, the solution appears to be converging to a configuration with two droplets per period 



(see Figure |T2D. All the solutions shown in this figure are numerically resolved; the small 
droplet is real rather than a numerical artifact, although it may later vanish as the solution 
evolves. 

• Random perturbations. 

Random zero-mean perturbations led to solutions with the same dynamics: perturbations 
of h = .5 relaxed to the mean, perturbations oi h = 1.05 relaxed to a translate of the positive 
periodic steady state h^s, and perturbations of h = 2 appeared to converge to a configuration 
of two droplets. 

4.5. q= 1.768. 

Characteristic features for q G (1.75, 1.79).' some positive periodic steady states are linearly 
stable, while others are linearly unstable; and there can be more than one positive periodic 
steady state with the same period and area. (See bifurcation diagrams |^-h, and [T^, §5.1].) 

We take n = 1, m = 1.768, q = 1.768, and compute solutions of 

4.5.1. q = 1.768. Perturbing the positive periodic steady states. For g-values in the interval 
(1.75, 1.794] (approx.) a new possibility arises: a heteroclinic connection between two fun- 
damentally different positive periodic steady states. We investigate this possibility in what 
follows. 

For Bond number B = 1.257 we consider two distinct positive periodic steady states, /issi 
and hss2, that have least period 27r, area A = 4.661, and have their local minima at x = 0. 
We denote the steady state that has lower minimum value by /issi, and the other by /iss2- 
Then we expect hgsi to be linearly stable and hss2 to be unstable, by [|19], Theorem 9] and its 



accompanying remarks, with hssi having lower energy. That is, hssi lies on the stable branch 
of the bifurcation diagram]^ and /iss2 lies on the unstable branch (since BP^l~'^A1~^ = 39.46.) 
Note also that the constant steady state h := hssi is linearly stable since Bh'^ ^ < 1. 

We consider even perturbations {£h"s)y odd perturbations {eh'^^), and random perturba- 
tions. Our simulations show that hssi has stability properties like the g = 1.5 steady state 
hss examined in § [4.4| . The other steady state /iss2 is unstable. For example, the initial data 
^ss2 + 10~^/igg2 yields a solution that converges to the constant steady state, as shown in 



Figure |T3|(a). The initial data /iss2 — ^^~'^Ks2 yields a solution converging to /igsi as t ^ oo, 
as shown in Figure 0(b). The observed behavior is very robust, and strongly suggests exis- 
tence of a heteroclinic connection from the unstable steady state hss2 to the stable one /igsi- 

Perturbing /iss2 with random zero-mean perturbations yields evidence of heteroclinic con- 
nections from hss2 to translates of /igsi and to the constant steady state. 
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Fig. 13: q = 1.768 and n 
/iss2 + 10^^/i"s2; solid: hss2- Solution relaxes to 
constant, (b) Dashed: /iss2 — 10~'*/iss2; solid: 
ft-ssi • Solution relaxes to hssi ■ 



1. (a) Dashed: Fig. 14: q = 2.5 and n = 1. (a) Dashed: hsa + 
10^"'/i"; solid: ft-as- Solution relaxes to ft-ss • (b) 



Dashed: - lO^^/i^;, 
in finite time. 



Solution touches down 



We now explain how we found the two steady states h^si and hss2 having the same period 



and area, since this is not completely obvious. First we plot E[a) for a G (0, 1), as in ||19 
Figure 5]. We seek ai < such that E{a\) = E{a2) with E'{ai) < and E'{a2) > 0. 
Choosing ai = .05, it is graphically clear from [|l^. Figure 5] that the desired a2 exists. To 
determine it, we compute ka^ and its period P(ai) = 6.703 and area A{ai) = 5.660, and 
hence E{ai) = 39.46. Taking a = ai and D = Di = 1 and Pss = 27t in the rescaling relations 
||T9| , eq. (24)] we evaluate the Bond number as i3 = 1.257. We then find six values of a such 
that E{a) < E{ai) at the first three values and E{a) > E{ai) at the last three values, so 
that a2 is between the third and fourth values. We interpolate E at these six values of a 
with a quintic polynomial and use Newton-Raphson iteration to find a2 = 0.5069 satisfying 
E{ai) = E{a2). We find has period ^(0:2) = 6.388 and area A{a2) = 6.115. Then 
P{a2) and the Bond number B determine the integration constant D2 = 0.8010 from |]T9|, eq. 
(24)], using again Pgs = 2n. Using Di and B and |1^, eq. (7)], we rescale ka^ to find /igsi; by 
similarly using D2 we rescale /cqj to find /iss2- Then we use /igsi and hss2 to determine finite 
difference steady states, as in 



4.6. q = 2.5. 

Characteristic features of q & (1.8,3).' positive periodic steady states are linearly unstable. 
'Mountain pass' scenario can occur. (See bifurcation diagram 0i and remarks after 119 , 
Theorem 11].) 

We take n = 1, m = 2.5, q = 2.5 and compute solutions of 

4.6.1. q = 2.5. Perturbing the positive periodic steady state. For Bond number B = 1.561 
we consider a positive periodic steady state hss with period 27r and area A^s = 4.335. This 
arises from rescaling ka with minimum height a = 0.2145 and period P^ = 6.869. Here hgs 
is linearly unstable. 
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Figure |l5(a) presents the solution with initial data hss + lO^^h"^, which converges to the 
constant steady state as time passes. Figure 0(b) shows the solution with initial data 
— 10~^^ss' which appears to touch down in finite time. These numerical results are 
qualitatively the same as for q = 0.5, in §4.2| . We found this behavior was very robust — 
we considered random perturbations (see § [4.1.1|) and found for each one that the solutions 
would either converge to the constant steady state or else touch down in finite time. 

4.6.2. q = 2.5. Perturbing the constant steady state. Just as in the q = —3 case, our numerics 
robustly confirm our predictions: for Bh'^ ^ < 1, small zero-mean perturbations of h yield 
solutions that relax to the constant steady state, while for Bh'^ ^ > 1 such perturbations 
yield solutions that appear to touch down in finite time. 

4.7. q = 4. 

Characteristic features of q & [3,oo).' positive periodic steady states are linearly unstable, 
and if a positive periodic steady state and a zero- angle droplet steady state have the same 



area, then the period of the former is less than the length of the latter. (See [jT8|, Theorem 7] 
and the proof of [ p!9| . Theorem 7].) 
We take n = l, m = 4, g = 4 and compute solutions of 

(9) ht = -{h^h^^^)^ - BQi^K 



''X Jx- 



This equation is 'super-critical' in the sense of Bertozzi and Pugh |^, since m > n + 2 [i.e. 
q > 3). According to the conjecture in [§, then, positive periodic solutions can blow up in 
finite time (||/i(-, t)||oo —* oo). Bertozzi and Pugh made the same conjecture for compactly 
supported weak solutions on the line, and proved blow-up can occur in finite time when 
n = 1 and m>n + 2 = 3 0. Specifically, they proved for such cases that if the compactly 
supported initial data ho{x) has negative energy 



dx < 0, 



then the compactly supported weak solution blows up in finite time, with its L°° and 
norms both going to infinity. 

Here we present computational evidence that smooth periodic solutions of (|^) can also 
blow up in finite time (there is no proof of this). Further, we find initial data that has 
positive energy yet still appears to yield finite-time blow-up, suggesting that negativity of 
the energy is not necessary for blow-up, in the periodic case. 

4.7.1. q = 4:, perturbing the positive periodic steady state. For Bond number B = 4.281 we 
consider a positive periodic steady state hss with period 2tt and area Ass = 3.213. This 
arises from rescaling ka with minimum height a = 0.2145 and period = 7.536. Here hss 
is linearly unstable and the energy £{hQ) = .01445 is positive. 
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• Even perturbations. 

We considered initial data hgs ± .OOl/igg. The initial data hgs + .OOlh"^ yielded a solution 
that relaxed to the constant solution as time evolved. The local extrema were fixed in space 
and, after a short transient, relaxed monotonically to the mean. 




Fig. 15: q = A,n= 1. Dashed: initial data hss — .001/i"g. In the figure, the maximum increases 
with time. The final profile, at t = 10.0282848199749, has 4096 meshpoints. 



The initial data hss — -OOl/igg yielded a solution that appears to blow up in finite time 
(see Figure ^). The extrema are fixed in space and, after a short transient, hmax increases 
monotonically towards infinity (the figure shows hmax increasing by a factor of 18.7). After 
a short transient, hmin decreases monotonically to a positive value as the singular time 
approaches. 

A self-similarity ansatz similar to the case q = —3 suggests that h(x, t) ~ (tc — t)^^^'^H({x — 
7r)/(tc — t)'^^^'^) as blowup approaches at x = vr; here tc is the time of blowup and if is a 
positive function with H(r]) ~ C|?7|^^/^ for large rj. Our computations are consistent with 
the above ansatz. Again, if we make the ansatz then we can estimate the blowup time tc, 
since taking the ratio of the computed values of h{7T, t) at two late times ti and t2 gives 
the value of {t^ - ti)"^/7(^c - ^2)"^^^ from which tc can be determined. We find that tr is 
slightly larger than the final resolved time. 

Self-similar blow-up for super-critical exponents has also been found for ht = —{h^hxxx)x — 
B{h^h,), in 01. 

• Odd perturbations. 

The odd perturbation 10~'^/igg yielded an evolution qualitatively like that shown in Figure 
T5| . The solution appears to blow up in finite time, with the location of the maximum moving 
in time. 

• Random perturbations. 

Random perturbations led to the same type of dynamics as seen with ±/i"g perturbations. 



4.7.2. q = 4, perturbing the constant steady state. Our numerics robustly confirm our pre- 



— g— 1 

dictions: for Bh < 1, small zero- mean perturbations of h yield solutions that relax to the 
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constant steady state, while for Bh'' ^ > 1 such perturbations yield solutions that appear to 
blow up in finite time. 

5. The effect of changing the mobility exponents, n and m 

In this section, we vary n and m in the equation ht = —{h^hxxx)x — B{h"^hrc)x- We think 
of these exponents as mobility parameters, since they determine the diffusion coefficients of 
the fourth and second order terms in the equation. When we change n and m, we keep 
q = m — n + 1 fixed, which means the steady states of the evolution are unchanged by 
§|2|. This allows us to ask three natural questions about the effects of changing the mobility 
exponents: 

1. Can a heteroclinic orbit between steady states be broken, or is it merely perturbed? 

2. Can the type of a singularity be altered {e.g. from finite-time to infinite-time)? 

3. Can the number oi singularities be altered {e.g. from one to two per period)? 

The next three subsections address these questions. 

Note that while one often thinks of the initial and terminal points of a heteroclinic con- 
nection as being isolated equilibrium points, here our equilibria are not isolated, for two 
reasons. First, the translates of a steady state are themselves all steady states with the same 
area. (Of course this translational freedom disappears when considering even perturbations 
of even initial data.) And second, when we consider configurations of droplet (compactly 
supported) steady states that consist of several droplets, the droplets can be translated, 
shrunk or expanded, subject only to the requirements that the total area (volume) be fixed 
and that the droplets remain disjoint. 

5.1. Perturbing heteroclinic orbits. For q = 2.5, we computed a 27r-periodic steady 
state with Bond number B = 1.561 and area 4.335, by rescaling with a = 0.2145 
and Pa = 6.869. The steady state is linearly unstable (see bifurcation diagram ^ with 
BPgl~'^A1~^ = 35.32.) We take initial data hss + .OOlh'^^; we expect to find solutions that 
relax to the constant steady state, hss = 0.6899. We vary the mobility exponents, by taking 
n = 0, 1, 2, and 3 in turn, and determining m from q = 2.5 = m — n + 1. That is, we 
compute solutions for the four evolution equations, all with the same initial data. 

We find that all four solutions relax to the constant steady state; the apparent heteroclinic 



orbit is robust under this change in mobility. In the top left plot of Figure |16|, we plot /imm 
and hmax versus time for the four solutions. The larger the exponent n, the longer it takes 
for the solution to relax to the constant steady state. 

Since hss < 1, hmin{t) < 1 for all time and there is a time ti, dependent on n, such 
that hmax{t) > 1 for t < ti and hmax{t) < 1 for t > ti. Near the minimum point, h^ > 
h^ > h"^ > h^, suggesting that the larger n is, the slower the diffusion will be (near the 
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Fig. 16: Dotted n — 0, dot-dashed n = 1, dashed n = 2, sohd n = 3. Left: q = 2.5. Top left: 
hmax{t) and hmm{t) versus t. Bottom left: hmax{t — ti) versus t — ti, where hmaxiti) = 1. 
Right: q = 0.5. Top right: hmax{t) and hminit) versus t. Bottom right: hmin(t — ti) versus 
t~ii where hmin{ii) = 1- 



minimum). Similarly, so long as the maximum is larger than 1 we have < < h"^ < 
near the maximum, suggesting that the larger n is, the faster the diffusion will be (near the 
maximum). This conflict of timescales appears to be mediated through the conservation 
of mass. Since the mean of the solution is conserved, we find that the solution moves as 
slowly as its slowest part (which is around the minimum): thus the larger n is, the slower 
the diffusion. This is demonstrated in the upper left plot of Figure |16]: the n = 3 solution 
takes longer to relax than the n = solution. Beyond the time ti there is no conflict; both 
the minimum and maximum should relax more slowly as n increases. We demonstrate this 



by plotting h^ax versus t — ti in the bottom left plot of Figure |T6 



For g = .5 we compute a 27r-periodic steady state hss with Bond number B = 0.9817 
and area Ags = 7.165. This arises from rescaling fc„ with minimum height a = 0.2145 and 
period Pa = 6.168. The steady state is linearly unstable (see bifurcation diagram [pD with 
BPgl~'^A1~^ = 36.29.) As for q = 2.5, we take initial data hss + .OOl/igg and exponents n = 0, 
1, 2, and 3. 

Again, we find that all four solutions relax to the constant steady state hss = 1-140; the 
apparent heteroclinic orbit is robust under this change in mobility. In the top left plot of 
Figure 0, we plot hmm and hmax versus time for the four evolutions. We see that the larger 
the exponent n, the longer it takes for the solution to relax to the constant steady state. 

Since hss > 1, hmax{t) > 1 for all times and there is a time ti, dependent on n, such that 
hmin{t) < 1 for t < ti and hminit) > 1 ioi t > ti. By the same logic as before, for t < ti, the 
time-scales will be dominated by the dynamics of hmm- This is demonstrated in the upper 
right plot of Figure |l^: the n = 3 solution takes longer to relax than the n = solution. 
Beyond the time ti, however, both the minimum and maximum should relax more quickly 
as n increases. We demonstrate this by plotting hmin versus t — ti in the bottom right plot 



of Figure 16. There we see the speeds of relaxation reverse, as expected. 
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For q = 2.5, we verify as follows that the profile of the solution is not largely affected by the 
mobility function h". First, we find the four 'half-times': the times at which hmin{ti/2) = -5. 
(The use of .5 is essentially arbitrary.) We find that the four solutions at their half-times 
differ by only 0.1% in the norm. We find analogous results for q = .5. This suggests 
that using hmin{t) to set a time-scale is an effective way of closely correlating two points on 
two different heteroclinic orbits. 

5.2. Changing the type of singularities. The choice of mobility coefficients in equation 
(|I|), hf = —{h'"'h^rcx)x — ^{h'^hx)x, affects whether a positive solution can become zero some- 
where in finite time. For example, if 3.5 <n<m<n + 2 then it cannot: the solution stays 
positive for all time H, §4.2]. (Note that m < n + 2 means q < 3.) On the other hand, if 
m > n + 2 then it is possible that the solution could blow up: ||/i(-, t) — oo in finite 
time. But even then we know from the methods of that if 3.5 < n < m then the solution 
remains positive as long as it exists. 

Here, we seek the critical exponent nc{q) such that if n > ndq) then positive initial data 
yield positive solutions for all time, while if n < ndq) it is possible for a positive smooth 
solution to touch down in finite time (becoming then a nonnegative weak solution). From 
above, if g < 3 then ndq) < 3.5. 

For q = 1, Goldstein et al. [0, §4] presented simulations with n = 1 that suggest a finite- 
time singularity is possible if i3 > 1. Bertozzi and Pugh presented numerical simulations for 
q = 1 and n = 3 in which the solutions remain positive for all time and appear to converge 
to one droplet per period as t — oo g]. This suggests that 1 < nc(l) < 3. 

Here we consider two further g-values, q = 2.5 and q = .5. In each case we take initial 
data hss — .OOlh"^, with the same steady states /igs as in ^5.1j . 

We saw for q = 2.5 and n = 1, in § [4.6| , that solutions appeared to touch down in finite 
time, hence 1 < nc{q) < 3.5. To further approximate the critical ra-value, we performed 
simulations with n = 1.25, 1.5, 1.55, 1.6, 1.65, 1.6625, 1.675, 1.7, 1.75, 2, 2.25, 2.5, 2.75, 3, 
and 4. Our findings suggest 

1.65 < ne(2.5) < 1.6625. 

In the left plot of Figure |1^ we plot log^Q hminif) versus t for most of these exponents. If 
hmin is decreasing at an exponential rate then the graph will be linear at large times. If hmm 
is decreasing to zero in finite time with an algebraic rate then the graph will go to — oo at 
some finite time, dropping down with a vertical slope. From the plot, iin = 2 (the rightmost 
graph) then hmin decreases monotonically in time, eventually decreasing with an exponential 
rate. For n = 1.75, 1.7, 1.675, or 1.6625, hmm decreases, then increases, and then ultimately 
decreases with an exponential rate. The solutions with n < 1.6625 appear to be touching 
down in finite time. However the n = 1.6625 simulation gives a note of caution; it is possible 
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that the simulations with n < 1.6625 would run until hmin became quite small but would 
then increase and ultimately decrease exponentially. Note: all of the simulations were run 
until the 8192 meshpoint simulation lost resolution, except for the n = 1.6625 simulation 
which required 65,536 meshpoints to resolve the solution when hmin was at its smallest. 
We ran the simulations many decades beyond those shown to verify the exponential rate of 
decrease. 




Fig.l7: logio(/i™„(i)) versus t. Left: q = 2.5; n = 1, 1.25, 1.5, 1.55, 1.6, 1.65, 1.6625, 1.675, 
1.7, 1.75, 2. Right: g = .5; n = 1, 1.5, 1.8, 1.85, 1.9, 1.95, 2. 



For g = .5 we similarly considered n = 1.25, 1.4, 1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.8, 1.85, 
1.9, 1.95, 2, 2.25, 2.5, 2.75, 3, and 4. We see very similar phenomena to the q = 2.5 case. 



The right plot of Figure |1^ is the analogue of the left plot and suggests that 

1.8 < n,(.5) < 1.85. 

For n > 1.85, the solutions stayed positive for the length of the computation and hmin 
decreased exponentially in time. 

5.3. Splitting singularities. Our work in § p.l| suggests that heteroclinic orbits can be 
preserved under some changes of the mobility. On the other hand, qualitative features of 
a solution can change significantly when the mobility is changed. For example, in § |5.2| we 
demonstrated that the mobility can affect the regularity of the solutions; for sufficiently large 
n, solutions are classical for all time while for smaller n, solutions can become nonnegative 
weak solutions in finite time. In this section, we demonstrate another effect of changing the 
mobility: a solution that touches down at one point per period can change into one that 
touches down at two. 

5.3.1. q = 2.5. We first consider q = 2.5 and n > 0. Then positive periodic smooth solutions 
of (|l]) remain bounded in for as long as they exist: || /;,(■, t) < M < oo by p|. We 
expect that these solutions will converge to a steady state, as t oo. But the positive 
periodic steady state is linearly unstable, and so we expect solutions to converge either to 
the constant steady state or to configurations of steady droplets (zero or nonzero contact 
angle). 
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As before, we take initial data hss — .OOl/igg with the same hss as in § ^.1| . In § |5.2| , we found 
that if > 1.6625 ~ nc{2.5) then solutions appear to stay positive for all time, with hmin 
decreasing to zero exponentially slowly in time. 

We find for n = 1 that the solution touches down in finite time at one point per period, 
consistent with a long-time limit of one droplet per period (see left plot of Figure |1^). For 
n = 2 the solution appears to be positive at all times and to touch down at two points per 
period in the long-time limit (see right plot of Figure |T8|). This suggests a long-time limit 
of two steady droplets per period. But it is impossible to contain two zero contact angle 
steady droplets in an interval of length 27r, as we argue shortly. In fact, we find that the 
small 'proto-droplet' is actually draining, with its maximum decreasing to zero like 
A similar phenomenon was observed by Constantin et al. |T2| §111, IV] with n = 1 and B = 
(and with different boundary conditions), although their proto-droplet seemed to decay like 
1/t. In our case, we find that the draining rate depends on n. 





Fig. 18: q = 2.5. Left: n ~ 1; hrnin{t) occurs at a; = 0. At all times there appears to be one 
droplet per period. Right: n = 2; at late times, the global minima flank a; = 0, suggesting 
a long-time limit of two droplets per period. But in fact the smaller droplet appears to be 
vanishing as hmin 0. 



Our simulations suggest that a second critical exponent, nc{q), governs the number of 
touch-downs per period, at least for the even perturbations we are using. If n < ficiq) then 
there appears to be one touch-down per period, occurring at x = 0. If n > hc{q) then 
there appear to be two touch-downs per period, with the position of the local minimum 
moving in time and with the solution being non-symmetric about the local minimum. The 
singularity splits as n increases through hc{q)- Goldstein et al. |T3|, §4C] observed something 
similar for ht = —{hhxxx)x ~ B{hhx)x [q = n = m = 1). Specifically, they found a single 
symmetric singularity that splits into a pair of asymmetric singularities as B increases from 
1 past B ^ 1.35. 

We find that 



1 < nc(2.5) < 1.25. 



32 



In the left plot of Figure |19] we plot the late-time profiles for a range of n. All five profiles 
shown are the final resolved solution with 8192 meshpoints. In the top plot, we plot the 
profiles from n = 1 and 1.25. The n = 1 profile has only one local minimum, while the 
n = 1.25 profile has two. In the bottom plot, we plot the profiles from n = 1.5, 1.6 and 
n = 1.7. Each profile has two local minima, with the distance between the minima increasing 
with n. 
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Fig. 19: Late-time profiles of h. Left: q = 2.5. Top: n = 1 (single minimum) and 1.25. 
Bottom: n = 1.5, 1.6, 1.7; width of proto-droplct increases with n. Right: q = .5. Top: 
n = 1 (single minimum) and 1.25. Bottom: n = 1.4, 1.6, 1.8, 2; width of proto-droplet 
increases with n. 



The n = 3 and 4 evolutions appear to be very similar to the n = 2 evolution. We did not 
do any simulations beyond n = 4 since the larger the value of n, the longer the simulation 
had to run before we could observe anything tangible — the mobility coefficient h"' is very 
small near the local minimum, where h is small. We could not compensate by taking large 
time-steps, because /i" could be quite large near the local maximum, since h^axit) > 1- 
(Large differences in time-scales are difficult to handle numerically.) 

We now prove our earlier claim that for q = 2.5 and B = 1.561, there cannot be two 
disjoint zero contact angle steady droplets in an interval of length 27r, if the total area of 
the droplets is A^g = 4.335 (the area of the initial data in § p.l| ). There does exist a single 
zero angle droplet steady state with that area and with length less than 2tt (the length is 
Pss = 4~s^ {Eo{2.5)/Bf = 5.287 < 2n), by ^ Theorem 7]. Hence there zs a zero contact 
angle droplet to which the solution on the right of Figure |l^ could relax. But if there were a 
pair of zero contact angle steady state droplets, with areas Ai = XAgs and A2 = {I — X)Ass, 
then the combined length of the two droplets would be 

Pi + P2 = (A-^ + (1 - A)-3) A;^ {Eo{2.5)/Bf . 

The righthand side is a convex function of A, achieving its minimum value 84.5905 at A = 1/2. 
This minimum value is greater than 271, and so one cannot fit the two droplet steady states 
in an interval of length 27r. Thus the simulations described earlier cannot be converging to a 
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pair of steady zero contact angle droplets. Indeed, our computations show the proto-droplet 
is slowly draining. 



5.3.2. q = .5. For g = .5 we took initial data hss — .OOl/igg (for hgs as in § |5.1| ) and observed 
phenomena very similar to those in the q = 2.5 case. We find there is an exponent nc(.5) 
such that if n < nc(-5) then there is only one touch-down per period and if ri > fic(-5) then 
there can be two. Again, our simulations suggest 

1 < nc(.5) < 1.25; 

see the plots to the right of Figure ITyi. In the right top plot, we present the final resolved 



solutions for n = 1 and 1.25 with 32,768 meshpoints. The local minima for n = 1.25 are 
much closer to each other than in the q = 2.5 plot. In fact, at 8192 meshpoints it appeared 
that there would be only one local minimum; as more time passed it split into two. In 



the right bottom plot of Figure |19| we plot solutions with n = 1.4, 1.6, 1.8, and 2. The 
n = 1.4, 1.6, and 1.8 solutions are with 8192 meshpoints, and the n = 2 solution is with 
16,384 meshpoints. The profiles are all resolved and were chosen to have comparable /imm- 
The distance between the local minima in the plot increase monotonically with n. 



The evolutions for n = 1 and 2 are very similar to those shown in the plots of Figure |18 
Specifically, for n = 2 the long-time limit appears to be one droplet. This is interesting since 
for g = .5 (unlike for q = 2.5), it is possible to have two disjoint zero contact angle steady 
solutions in an interval of length 2tt, as we now show. We have a Bond number B = 0.9817 
and area A^s = 7.165 (from §p.l|). A single zero contact angle steady state would satisfy [jl^ 
§3.1] 

where Eq{.5) = 32.86. We find Pgs = 6.038 < 27r. Thus a single zero contact angle steady 
state is a potential long-time limit. For two zero contact angle steady states, we find 

P, + P,= (A^/^ + (1 - Xy/') iEoi.5)/Bf' . 

This is a concave function with its maximum at A = 1/2. We find that Pi + P2 < 2tt if 
A < 10^"^ (approx.) and so one can have two zero contact angle droplets — but one of 
them must be fairly small. For example, if A = 10~^ then the length of the smaller droplet 



is Pi = 0.2404. In the bottom right plot of Figure |T9| the profile for n = 2 (the longest 
proto-droplet shown) has length .2431, which is close to .2404; thus a two-droplet long-time 
limit is at least a possibility. However, like for q = 2.5, the proto-droplet appears to shrink 
in time, with its maximum value decreasing to zero like (As for q = 2.5, the draining 

rate also depends on n.) We did not succeed in finding a q, n, and initial data that yield 
a solution with a multi-droplet configuration as its longtime limit. But we believe it should 
be possible to find this somehow. 
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6. Numerical methods 

The numerical simulations are done using a finite-difference evolution code. Throughout 
this section, the diffusion coefficients are represented as functions / and g\ we use power law 
coefficients f{y) = and g{y) = By^ in our simulations for this paper. The exponent ^ 
represents the t"^ time-step and is the numerical approximation of the solution at time 
t(, = iAt. 

6.1. The evolution code. We use an adaptive time-stepping scheme based on a Crank- 
Nicolson scheme: 

Here the diffusion coefficients are evaluated at /i^^^/^, which we find by linearly extrapolating 
the solutions h^~^ and to time ti + At/2. The scheme is O(At^) on time intervals where 
At is fixed and is 0{At) at the time when the timestep is changed. We explain the adaptive 
timestepping in § |0[ 

Finding /i^"*"^ reduces to solving a linear problem, which we write in a residual formulation: 
(10) Cz = -At [/(/.^^^/^)(/.L. + 

where z = h^^^ — and r{y) = g{y)/ f{y), and the linear operator £ is defined by 

Cz:=z + At^ [f{h'^'/'){z,.. + r{h'-''/')z,)]^. 
The scheme uses h^~^ and to compute h^~^^; for the first step we take = h!^. 



We now perform a linear stability analysis of the scheme about the constant steady state 
h = h, for power law coefficients f{y) = and g{y) = Ey"^. The linearized equation is 
ht = —h'^hxxxx ~ Eh^hxx- Initial data h^[x) = ecos{kx) yields /i^^^ = a{kY^^ ecos{kx) 
where 

i_^7rp(fc2_g7r"-") _ i-^(fc) 
^ " 1 + f 7rp(p _ gTT-) ~'- 1 + fi{k) • 

Hence there is a band of unstable modes: if < /c^ < Bh"^ " {i.e. fi{k) < 0) then |cr(A;)| > 1 
and the initial data hQ{x) = ecos{kx) yields a solution that grows exponentially in time. 
We consider a numerical scheme linearly stable if perturbations outside this band are not 
amplified: 

k^ > kl ■= Bir~" =^ \a{k)\ < 1. 

It follows immediately from the form of the growth factor that the Crank- Nicolson scheme is 
linearly stable, as expected. Such a linear stability analysis provides a useful guide, though 



it is directly relevant only for small perturbations of fiat steady states. 
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For the spatial discretization, the key issue is to implement the scheme in a way that 
preserves steady states. By (H), a steady state hss satisfies hssxxx + r{hss)hssx = 0. Such an 
'analytic steady state' will not generally be a 'finite-difference steady state', although hgs will 
be 0(Aa;^)-close to the finite-difference steady state h satisfying the following discretization: 

(11) - + 3k - hi-i + Ax=^^^^^^^:i^^tlM (h^_^^ -k)=0, i = l...N. 

The meshpoints are xi = Ax, X2 = 2Ax, ■ ■ ■ , xn = X, where X is the length of the interval, 
and we denote the function values at the meshpoints with subscripts: hi, h2, . . . Hn- The 
function is periodic: Hq = h^. 

To implement the residual formulation (p!0|), we apply the OlAx"^) approximation 

fi- I ^i+l ~ + 3Zj_i — Zj_2 Zi — Zi-i\ 

^ ' ~A^[ A^^ + Ax J ' 

where the subscripts i+ and i— denote the right-average and left-average, for example: 

ri+.- 2 ' 2 

This approximation yields a 0(Ax^)-accurate matrix approximation C of the operator C. 
Using this, the residual formulation is written 

Cz = RHS{h^'\h^). 

The N X N matrix C is pentadiagonal periodic. The righthand side of (^) is discretized 
analogously. 

It follows immediately from the condition (^) for a finite difference steady state that 
the above time-stepping scheme preserves finite-difference steady states. To ensure this we 
factored / out in (|T^) before doing the finite-difference approximation of f{h)hxxx + g{h)hx, 
because the relation fr = g does not hold in the discrete setting: fi+Ti^ ^ gi+ for example. 
Also, by factoring / out we have also isolated the pressure gradient term, {hxx + B/(lh'^)x, in 
the equation. 

6.2. Timestepping and accuracy. The adaptive timestepping controls the accuracy as 
follows. An error tolerance is set, e = 10"^^. At each time step, we first use the Crank- 
Nicolson scheme to compute hi, an approximation of the solution at time t + At. We then 
take two timesteps with At/2 to compute /12, another approximation of the solution at time 
t + At. For some constant C, the error is bounded [0, §5.2] by the difference of hi and h2: 

\\h{-,t + At) -/islloo < - /12II00. 
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li \\hi — h2\\oo > £ then we replace At with At/2 and try again (without advancing in time). 
If \\hi — h2\\oo < b/10 then we replace At with 2At and try again. If \\hi — /i2||oo lies between 
e and e/10 then we take the solution at time t + At to be /i2. 

Adaptive timestepping takes at least three times longer than using the Crank-Nicolson 
scheme with a fixed timestep. For this reason, we performed all of our exploratory studies 
using a fixed time-step. Once we found phenomena of interest, we re-ran using adaptive time- 
stepping. Since the first time-step has O(At^) local truncation error, rather than 0(At'^), 
the adaptive time-stepper initially refines At to meet the tolerance. Also, most runs had an 
initial fast transient (see § |4.1.1| ) which required early refinement of the timestep. 

Admittedly, since we do not know the constant C this error control is valid only as long 
as C is not large. In practice we find that, after the initial transient, the timestep is rarely 
reduced, except near times when the run has to be stopped anyway in order to point-double. 



6.3. Computing the finite-difference steady state. Here, we only discuss the case of 
power law coefficients. Given uniformly distributed meshpoints between and X = 27r, 
we seek a finite-difference steady state h that solves the A^ equations (p!l| ) to the level of 
round-off error. 

We solve the A^ equations ([Tl| ) simultaneously using Newton-Raphson iteration. To do 
this, we need a good first guess for h. In the following, we describe how we find a first guess 
and then how we execute the Newton-Raphson iteration. 

Given the exponent g 7^ we first compute the rescaled steady state k = at the A^ 



points X = P/N, 2P/N, ■ ■ ■ ,P. As described in ||T7| , §6.1], we do this by viewing the steady 
state equation 

(13) fc.a= + ^^ = 0, k{0) = a, fc.(0)=0, 

as an initial value problem in x. We verify that k is spectrally accurate by using a discrete 
fast Fourier transform to check that the power-spectrum is fully resolved. Since for numerical 
purposes, k is an exact solution of equation (0), in the following we refer to it as an 'analytic 
steady state'. 

Once the analytic steady state k is known, we rescale it to find an analytic steady state 
hss of period 27r, by taking £) = 1 in the rescaling eq. (7)] and defining 



B = { — , and hss = { — \ k. 

This gives hss at A^ meshpoints. By construction, hss is an analytic steady state for {hss)xx + 
{Bhl - 1)/? = with period 27r. 
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To find the finite-difference steady state h close to /igs we need to solve the equations 
(0), which we write as 

F{h) ■= Mh + V{h) = 

where M is a tetradiagonal periodic matrix and V{h)i is a nonlinear function of hi and /ij+i. 
The Newton-Raphson iteration is 

j^new ^ j^oid _ {^DF(h"^'^))-^F{h°^'^). 

To iterate, one has to solve DFx = F{h"^'^). We find that DF is a singular matrix of rank 
— 1. We solve DFx = F(h°^'^) using the singular value decomposition of DF obtained 
using LAPACK's 'dgesvd.f ' to solve for x. The iteration is then started with the initial guess 
hss and stopped when the largest error in the A^ equations (|TT]) is less than 10"^''. 

In a different approach, one might attempt to compute the finite difference steady state 
h with a relaxation method, by computing the finite-difference solution of the evolution 
equation ht = Ah + {Bh^ — \)/q for large t. This seems unlikely to succeed because steady 
states of this PDE have the same linear stability properties (with respect to zero-mean 
perturbations) as the steady states we are trying to compute, making it extremely difficult 
to obtain convergence to linearly unstable steady states. 

Note: For g = 1, the finite-difference steady states satisfy a linear problem: 

hi+2 — 3/ii+i + 3hi — hi_i + Ax^B ^/ij+i — h^ =0, % = \ . . .N . 



In §0, we perturb a nontrivial 27r-periodic finite-difference steady state. This will be close to 
an analytic steady state. Even analytic steady states are a + 6cos(v^a^) and are 27r-periodic 
if ^JB is an integer. Sampling such a steady state on a uniform mesh gives 

hj = a + b cos (VBjAx), 

and one can check that h is a finite-difference steady state provided 

B = 2 '-Z^f^"^ =B-OiAx^ 
Ax^ 

That is, there are nontrivial analytic steady states for a countable collection of Bond numbers, 
i3 G {1^, 2^, 3^, . . . }, and nontrivial finite difference steady states for a nearby countable set 
of Bond numbers B. 

6.4. Stopping criteria and issues for singularities. To test whether to stop the code, 
we compute the minimum value of h^ at each time-step. If this minimum is ever less than 
or equal to zero, we stop the code. As discussed in §§^ and |, we find that the stopping 
criterion is indeed often met. While this suggests a finite-time singularity, we emphasize 
that the code has not been written in a way to carefully resolve such singularities — the 
singularity may occur in infinite time, with the stopping criterion being met in finite time 
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because of instabilities causing oscillations in the profile. In practice, the stopping criterion 
was always met after the solution became spectrally unresolved. 

The code was designed to preserve the periodic steady states, rather than being written to 
preserve positivity. (There are a number of different approaches to ensuring positivity: we 
refer interested readers to JI], |2^.) For this reason, it might be that the code is stopping 
spuriously. Also, as the code has no local mesh refinement, we have to over-resolve much of 
the solution in order to resolve the solution where it is tending to zero. This over-resolution 
away from the singular points slows the computation significantly. 

That being said, we use the code primarily to study nonsingular behavior. Some of our 
results are suggestive of finite-time singularities; we present them with the above caveats. To 
find fine details of the temporal and spatial scales of the singularities, we would implement 
a code that preserves positivity and has an adaptive spatial mesh. 

Finally, we are not using a fully-implicit timestepping scheme; a basic such scheme would 
only be 0{At) but could be made O(At^) using Richardson extrapolation [T^. It would 
likely be very stable, but slow. For speed, we chose a partially implicit scheme and then 
checked for numerical instabilities when we held the timestep fixed. We observed none. 
Also, it seems unlikely that our adaptive timestepper was taking small timesteps in order 
to control numerical instabilities, since after the initial transient, the timesteps were refined 
only when the solution was becoming singular. 



7. Conclusions and Future Directions 

We have shown numerically for the evolution equation ht = —{h'"'hxxx)x — B{h^hx)x that 
our linear stability theorems in [|I8| accurately predict the short-time nonlinear behavior 
of the solutions near positive periodic and constant steady states. We have found strong 
evidence for the existence of heteroclinic connections between steady states, as suggested by 
our theorems on the energy levels of steady states [|T9] . We have further observed a mountain 
pass scenario in which perturbations of a periodic steady state relax towards either a droplet 
or a constant steady state. 

All of this suggests that the energy landscape through which the solutions travel is fairly 
simple, and that understanding the relative energy levels of the steady states gives consid- 
erable insight into that landscape. 

It is worth recalling that the evolution equation (P describes gradient flow for the energy 
£ defined in (|^), with respect to the following weighted inner product. Let h{x,t) be a 
positive smooth function that is X-periodic in x, and for each t, define an inner product on 
functions a, 6 G C{Tx) having mean value zero by 

{a,b) := [ A'{x)B'{x)h{x,tYdx, 
Jo 
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where A,B e C^{Tx) satisfy = a and (h'^B^)^ = b. Then the equation (g) is 

equivalent to: 

SS 

— = ~{ht, (p) for all G C{Tx) having mean value zero. 

Hence the variation of the energy is most negative in the direction (f) = hf, so that the 
evolution equation for h{x, t) simply describes flow by steepest descent on the energy surface 
of S, with respect to the inner product (, ). Note that this inner product is time- dependent 
since it depends on /i", i.e. on the solution itself. 

The above gradient flow formulation with weighted inner product was observed by 
Taylor and Cahn ||2^; their evolution (7a) contains our equation (0). For the special case of 



the Cahn-Hilliard equation the inner product is unweighted, since the fourth order term in 
the equation is linear. Gradient flow ideas for related equations have been used in 0, ^ 



and a Wasserstein-fiow idea in 22 



In §^ we presented numerical results on the persistence of heteroclinic connections under 
changes in the mobility parameters n and m. There we changed n and m in a way that 
preserved q = m — n + 1, and thus preserved the energy S and also the steady states (which 
are critical points of the energy). Hence our change in mobilities does not change the energy 
landscape. But it does change the weight /i" appearing in the inner product (a, b) and in 
the equations for A and B, and this is how changing the mobility affects the evolution. In 
§^ we found the timescale of the solution changed noticeably in response to changes in the 
mobilities, even though the shape of the solution changed little. 

Lastly, in §|] we further investigated critical mobility exponents, such as the critical n 
above which solutions remain positive for all time (in other words, the critical exponent 
for film rupture or pinch-off). An interesting question for the future is to find formulas for 
the critical mobility exponents. These critical exponents determine important qualitative 
features of the evolution and determining them would shed considerable light not only on 
the equation (|l]) studied here, but also on related equations that arise from physical models. 



Acknowledgments. Laugesen was partially supported by NSF grant number DMS-9970228, 
and a grant from the University of Illinois Research Board. He is grateful for the hospitality 
of the Department of Mathematics at Washington University in St. Louis. 

Pugh was partially supported by NSF grant number DMS-9971392, by the MRSEC Pro- 
gram of the NSF under Award Number DMR-9808595, by the ASCI Flash Center at the 
University of Chicago under DOE contract B341495, and by an Alfred P. Sloan fellowship. 
The computations were done using a network of workstations paid for by an NSF SCREMS 
grant, DMS-9872029. Part of the research was conducted while enjoying the hospitality of 
the Mathematics Department and the James Franck Institute of the University of Chicago. 



40 



Pugh thanks Todd Dupont and Bastiaan Braams for illuminating conversations regarding 
numerical issues. 

References 

J. W. Barrett, J. F. Blowey, and H. Garcke. Finite element approximation of a fourth order nonlinear 

degenerate parabolic equation. Num,er Math, 80(4):525~556, 1998. 

E. Beretta, M. Bertsch, and R. Dal Passo. Nonnegative solutions of a fourth order nonlinear degenerate 
parabolic equation. Arch Ration Mech Anal, 129:175-200, 1995. 

F. Bcrnis and A. Friedman. Higher order nonlinear degenerate parabolic equations. J Diff Eq, 83:179- 
206, 1990. 

A. J. Bernoff, A. L. Bertozzi and T. P. Witelski. Axisymmetric surface diffusion: dynamics and stability 

of self-similar pinchoff. J Statist Phys, 93(3 4), 725 776, 1998. 

A. L. Bertozzi, M. P. Brenner, T. F. Dupont, and L. P. Kadanoff. Singularities and similarities in 
interface flow. In L. Sirovich, editor, Trends and Perspectives in Applied Mathematics, volume 100 of 

Applied Mathematical Sciences, pages 155-208. Springer- Verlag, New York, 1994. 

A. L. Bertozzi and M. Pugh. The lubrication approximation for thin viscous films: the moving contact 
line with a 'porous media' cut off of van der Waals interactions. Nonlinearity, 7:1535-1564, 1994. 
A. L. Bertozzi and M. Pugh. The lubrication approximation for thin viscous films: regularity and long 
time behavior of weak solutions. Com,mun Pur Appl Math, 49(2):85-123, 1996. 

A. L. Bertozzi and M. C. Pugh. Long wave instabilities and saturation in thin film equations. Commun 
Pur Appl Math, 51:625-661, 1998. 

A. L. Bertozzi and M. C. Pugh. Finite-time blow-up of solutions of some long-wave unstable thin film 
equations. Indiana Univ Math J, to appear, 2000. 

A. L. Bertozzi and M. C. Pugh. Presented at 1997 APS Division of Fluid Dynamics meeting. 
J. Bricmont, A. Kupiainen, and J. Taskinen. Stability of Cahn-Hilliard fronts. Commun Pur Appl Math, 
52(7):839-871, 1999. 

P. Constantin, T. F. Dupont, R. E. Goldstein, L. P. Kadanoff, M. J. ShcUcy and S.-M. Zhou. Droplet 
breakup in a model of the Hele-Shaw cell. Phys Rev E, 47(6):4169-4181, 1993. 

R. Goldstein, A. Pesci, and M. Shelley. Instabilities and singularities in Hele-Shaw flow. Phys Fluids, 
10(ll):2701-2723, 1998. 

G. Griin and M. Rumpf. Nonnegativity preserving convergent schemes for the thin film equation. Numer 
Math, to appear, 2000. 

G. Griin and M. Rumpf. Simulation of singularities and instabilities arising in thin film flow. Submitted 
to Euro J Appl Math, 2000. 

A. Iserles. A first course in the numerical analysis of differential equations. Cambridge University Press, 
Cambridge, 1996. 

R. S. Laugesen and M. C. Pugh. Properties of steady states for thin film equations. European J Appl 
Math, to appear, 2000. 

R. S. Laugesen and M. C. Pugh. Linear stability of steady states for thin film and Cahn-Hilliard type 
equations. Arch Ration Mech Anal, to appear, 2000. 

R. S. Laugesen and M. C. Pugh. Energy levels of steady states for thin film type equations. Preprint, 
2000. 

P. Manneville. Dissipative structures and weak turbulence. Academic Press Inc., Boston, MA, 1990. 
A. Oron and S. G. Bankoff. Dewetting of a heated surface by an evaporating liquid film under conjoin- 
ing/disjoining pressures. J Colloid Interface Sci, 218:152-166, 1999. 

F. Otto. Lubrication approximation with prescribed non-zero contact angle: an existence result. Com- 
mun PaH Diff Eq, 23:2077-2164, 1998. 

J. E. Taylor and J. W. Cahn. Linking anisotropic sharp and diffuse surface motion laws via gradient 
fiows. J Stat Phys, 77(1-2):183-197, 1994. 

M. B. Williams and S. H. Davis. Nonlinear theory of film rupture. J Colloid Interf Sci, 90(l):220-228, 
1982. 

T. P. Witelski and A. J. Bernoff. Stability of self-similar solutions for van der Waals driven thin film 
rupture. Phys Fluids, ll(9):2443-2445, 1999. 

T. P. Witelski and A. J. Bernoff. Dynamics of three-dimensional thin film rupture. Preprint, 2000. 



41 



[27] W. W. Zhang and J. R. Lister. Similarity solutions for van der Waals rupture of a thin filni on a solid 

substrate. Phys Fluids, ll(9):2454-2462, 1999. 
[28] L. Zhornitskaya and A. Bcrtozzi. Positivity preserving schemes for lubrication-type equations. SIAM J 

Numer Anal, 37(2):523-555, 2000. 

EMAIL CONTACT: laugesen@math.uiuc.edu, mpugh@math.upenn.edu 
Department of Mathematics, University of Illinois, Urbana, IL 61801 
Department of Mathematics, University of Pennsylvania, Philadelphia, PA 19104 



